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Abstract 


Two simulations, one batch and one real-time, of an 
acroelastically scaled wind-tunnel model were developed. 

The wind-tunnel model was a full-span, frec-to-roll model of 
an advanced fighter concept. The batch simulation was used 
to generate and verify the real-time simulation and to test 
candidate control laws prior to implementation. The real-time 
simulation supported hot-bench testing of a digital controller, 
which was developed to actively control the elastic 
deformation of the wind-tunnel model. Time scaling was 
required for hot-bench testing. The wind-tunnel model, the 
math models for the simulations, the techniques employed to 
reduce the hot-bench time- scale factor, and the verification 
procedures arc described. 
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ilh control mode due to klh time derivative of 
elastic mode 

[A],[B] general continuous system dynamic matrices 
[C],[D] 


l wing chord, 39.76 in 

F 1 ^ total generalized force on elastic (flexible) mode i 
[F], [G] general discrete system dynamic matrices 
[Go] discrete system input matrix applied to [u k ] 

[Gj] discrete system input matrix applied to (u k+1 } 

[‘G*J diagonal modal damping matrix, G f (i,i) = .03 
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h integration step size, seconds 

1 iMg control surface hinge moment, in-lbs 

k§ actuator steady state gain 

f f 2 

f K . ] diagonal modal stiffness matrix, K (i,i) = m^ 

|M f J diagonal modal mass matrix, M f (i,i) = m x 

[M cf ] matrix of control mode to elastic mode inertial 

coupling 

q dynamic pressure, lbs/in 2 

rl rate limit 

s Laplace variable 

{ u k } vector input at time t=kh 

V airspeed, in/sec (6696 in/sec at Mach = .5) 

x stale 

{ x k } state vector at time t=kh 

8[ generalized coordinate - control mode i 

damping of second-order denominator term of 
typical actuator transfer function, rad/scc 
Tjj generalized coordinate - elastic mode i 

output of Drydcn turbulence transfer function 
p density, 1.146xl0* 7 lbs-sec 2 /in 4 or slinches/in 3 

<jg RMS turbulence intensity, in/sec 

u Gaussian random number 

co$ frequency of second-order denominator term of 

typical actuator transfer function, rad/sec 
C 0 g break frequency of turbulence transfer function, 

(2tc)( 17.32) rad/sec 

01^ break frequency of anti-aliasing transfer function, 
(2 k)( 25.0) rad/scc 
[] matrix 

( } column vector 


Subscripts 

a aerodynamic 

aa anti-aliasing 

as anti -symmetric 

c commanded 

g quantity associated with turbulence 

lag aerodynamic "lag" quantity 

lin linear 

neg negative 

pos positive 

ni no load 

STALL actuator dynamic stall value 
sy symmetric 

0 quantity associated with position 

1 quantity associated with rate 

2 quantity associated with acceleration 

8 quantity associated with actuator 
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S u pers cri pts 

f associated with elastic (flex) mode equations 

ff effect on elastic due to elastic (flex-flex) 

fc effect on elastic due to control (flex-control) 

fg effect on elastic due to gust input (flex -gust) 

c associated with control mode equations 

cf effect on hinge moment due to flex (control -flex) 
cc effect on hinge moment due to control 

eg effect on hinge moment due to gust input 

Abbreviations 

AB2 Adams-Bashforth second-order predictor 

AFW Active Flexible Wing 

CAMAC Computer Automated Measurement and Control 
ISAC Interaction of Structures, Aerodynamics, and 

Controls 

LaRC Langley Research Center 

LEI Leading Edge Inboard 

LEO Leading Edge Outboard 

PPU Peripheral Processor Unit 

psf pounds per square foot 

RK2 Runge-Kutta second-order predictor-corrector 

RMS root-mean-square 

RVDT Rotary Variable Differential Transducer 

TOT LaRC Transonic Dynamics Tunnel 

TEI Trailing Edge Inboard 

TEO Trailing Edge Outboard 

Introduction 

The simulations described in this paper were developed as 
part of the ongoing Active Flexible Wing Wind-Tunnel Test 
Program ,1L231 a collaborative effort between NASA Langley 
Research Center and Rockwell International Corporation. 
Three wind-tunnel tests have been completed and the final test 
is scheduled for February 1991. The program objective is the 
validation of analysis and synthesis methodologies as applied 
to the multi-function control of a sophisticated aeroelastic 
wind-tunnel model The control functions being investigated 
include suppression of flutter, roll performance maximization, 
and load alleviation. Only the simulation models developed to 
support flutter suppression are discussed in this paper. 

Flutter is a potentially explosive dynamic instability that can 
occur when a sufficiently flexible wing begins to extract 
energy from the fluid stream. During the most recent tunnel 
entry, completed in November 1989, flutter suppression was 
successfully demonstrated at a dynamic pressure 24 percent 
above a measured open-loop flutter point. P* 3 ! 

Simulation Roles 

Two distinct, but interrelated, simulations were developed 
to support the AFW test program, a batch simulation and a 
real-time simulation. The batch simulation served as a "truth" 
model, and was used to: (1) evaluate the control laws to 
predict performance and establish gain and phase margins; 

(2) provide models and data files for the real-time hot-bench 
simulation; and (3) verify the real-time simulation. 

The real-time simulation of the model/wind-tunnel 
environment served to verify the functionality of the digital 
controller system in a hot-bench laboratory. End-to-end 
verification and debugging of the complex, one-of-a-kind 
digital controller system was essential, since whenever flutter 


testing is undertaken, both the model and the tunnel are at 
risk. 


Wind-Tunnel Model 

Figure 1 shows the AFW model mounted in the LaRC 
Transonic Dynamics Tunnel during the November 1989 test. 
The sting mount has an internal ball bearing arrangement that 
allows the model to roll +/- 145° about the sting axis. The 
fuselage is connected to the sting with a hydraulically 
powered pivot so that the model can be remotely pitched from 
approximately -1.5° to +13.5°. Destabilizing mass ballast 
was added to each wingtip so that the model would flutter 
within the operating envelope of the TOTf^l The tip ballast 
serves to lower both the first-bending and the first-torsion 
elastic mode frequencies, with the predominant effect on the 
first-torsion mode. The result is that the first-torsion and the 
first-bending elastic modes combine to form the primary 
flutter mechanism at a lower dynamic pressure than was the 
case for the original model (with no tip ballast). The tip 
ballast can also be rapidly decoupled in pitch from the wingtip 
by releasing a hydraulic brake. When decoupled, the tip 
ballast is restrained in pitch by a soft spring. Decoupling the 
tip ballast proved to be an effective flutter stopper during 
testing, providing a significant safety margin. For the flutter 



Figure 1 .- AFW model mounted in the LaRC Transonic 
Dynamics Tunnel. 



Figure 2.- Instrumentation of the AFW model. 
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suppression tests conducted in the November 1989 test, the 
roll brake was engaged, so that the model was not free to roll. 

Figure 2, a drawing of the model, illustrates the locations 
of the control surfaces and selected instrumentation. There 
are 8 control surfaces and 13 accelerometers. With the 
addition of the tip ballast, two wing accelerometers, 
previously co-located with the leading edge inboard control 
surfaces, were moved to the tip ballasts. Reference 4 
describes in detail the AFW wind-tunnel model prior to 
adding the tip ballast. The input/output signal list for both the 
wind-tunnel model and the supporting simulations is given in 
Appendix A. 

Hot-Bench Laboratory 

The AFW hot-bench simulation set-up, depicted 
schematically in figure 3, utilized the central real-time facility 
at LaRCt 5 !. The LaRC real-time facility consists of nodes or 
simulation sites that communicate by means of a 50-mcgabii- 
per-second fiber-optic digital-data network. The various 
simulation nodes on the network included two Control Data 
Cyber 175 computers, engineering control consoles, various 
aircraft cockpits, and motion base hardware. For the AFW 
hot-bench simulation, one Cyber 175 was used to integrate 
the equations of motion. An Adage graphics computer, used 
in this study, communicates directly with a Cyber 175 
through a port. New Terabit Eagle 1000 graphics computers, 
currently being installed at LaRC, will communicate over the 
fiber-optic network as another simulation node. Communi- 
cations with the digital controller occurred over analog lines in 
the same manner as when the controller was connected to the 
wind-tunnel model at the LaRC TDT. 

Real-Time Graphics 

To assist in visualizing the simulated model dynamic 


behavior, a real-time display (figure 4) was developed on an 
ADAGE graphics computer. The display presents model 
pitch and roll, control surface deflections, and total structural 
deformation of the simulated wind-tunnel model. The display 
is based on a finite-clement structural model of the test article. 
A subset of the finite-element structural nodes were used to 
highlight the main body, wings, and the eight control 
surfaces. The dashed line represents the undeformed 
configuration. The deformations can be exaggerated for ease 
of viewing by a factor set at the console. 

Simulation Time Scaling 

The AFW hot-bench simulation operated at a time scale 
slower than 1:1 (real time). Time-scale is a function of the 
integration step (h) and the computer frame (T). If T is larger 
than h, the simulation runs at l:(T/h) "slow." Since there was 
no human operator in the hot-bench loop, a time scale other 
than 1:1 could be accommodated. Several factors prevented 
the AFW hot-bench simulation from operating in real time. 
The control computer was scheduled to run at 200 Hz in the 
wind-tunnel. To avoid excessive digitally induced time 
delays, the hot-bench simulation needed to update at least 
twice for each digital controller frame, requiring a 400 Hz rate 
for the simulation if it were to run in real time. The minimum 
frame lime available on the Cyber 175 real-time clock was 5 
milliseconds (200 Hz). The simulation model itself was 
sufficiently complex that a computational frame of at least 
12.5 milliseconds (80 Hz) was required. Furthermore, since 
there are only two Cyber real-time computers and many real- 
time jobs, the hot-bench simulation often shared a Cyber 175 
with another job. With only one half of the Cyber computing 
power available, the 80 Hz simulation update rate would be 
further reduced to 40 Hz. Time-scale must, therefore, be an 
easily adjusted parameter for any dynamic component in the 
hot-bench loop. 
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Figure 4.- Adage generated real-time wireframe image of the simulated AFW model. 


Digital Controller Apparatus 

The control computer consisted of a SUN 3/160 
workstation augmented with a digital signal processor, a 
floating point array processor, and analog/digital conversion 
boards.! 2 * 6 ^ During both hot-bench and tunnel testing, signals 
to and from the digital controller go through a Rockwell- 
NASA interface box. The interface box has several functions, 
one of which is to house anti-aliasing filters that condition the 
signals from the wind-tunnel model before being sampled by 
the digital controller. Because of the time scaling, the anti- 
aliasing filters in the interface box were bypassed during hot- 
bench testing (figure 1). The anti-aliasing filter dynamics 
were included in the hot-bench simulation. Since the control 
laws were digitized assuming an update rate of 200 Hz in the 
wind tunnel, when the hot-bench simulation ran at 1 :n slow, 
the controller was clocked at (200/n) Hz to be dynamically 
equivalent. 

In addition to its primary function of implementing a 
selected control law at 200 Hz, the digital controller 
performed a variety of support functions. Static checks, 
dynamic data acquisition and storage, and sine-sweeps were 
performed. Data from 20 second sine-sweeps were shipped 
to another computer wherein both open-loop and closed-loop 
plant estimations were performed.! 7 ! 

Simulation Math Model 

Linear aeroelastic descriptions for the symmetric and anti- 
symmetric boundary conditions were generated using a set of 
aeroservoelastic design and analysis programs developed at 
the LaRC called ISAC, short for "Interaction of Structures, 
Aerodynamics, and Controls"! 8 !. Doublet lattice aerodynamic 


theory was used. These models were combined with 
empirical data to form a "whole aircraft" model wherein left 
and right actuators were modeled individually. 

Actuators 

Frequency responses for the eight individual actuators 
were measured with the wing elastic motion restrained. In the 
frequency range of interest, third-order transfer functions, 
with parameters optimized in a least square sense, produced 
good fits with measured frequency response data. In general, 
right and left members of an actuator pair required different 
parameters to achieve a good fit and were so modeled. 

All the actuator transfer functions had the following form, 

k s ag co 2 , 

S(s) = l l 8 

6 c (s) (s+a 8 )(s 2 + 2£ 8 co 8 s + wj) 

where kg is the steady state gain, ag is the first-order pole 
location, £g is the damping of the complex pair, and cog is the 
frequency of the complex pair. The second-order complex 
pair results from the compressibility of the hydraulic fluid 
together with compliance of the structure. The first-order pole 
reflects the dynamics of hydraulic fluid flowing through a 
small orifice whose size is regulated by position feedback 
error. Rate limits, as a function of load, were specified by the 
manufacturer and result from the maximum rate that hydraulic 
fluid can pass through a small orifice for a given difference 
between supply and chamber pressures. The first-order pole 
part of the transfer function was then a reasonable place to 
apply rate limits. An initial, linear rate was first calculated as 
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Figure 5.- Turbulence intensities used for simulation. 


follows, 

*°/w = a P ( 5 c' Xo) 


surfaces were obtained they were resolved into symmetric and 
anti-symmetric components that became inputs into the 
acroclastic equations. 


Positive and negative rate limits based on no-load rate limits 
modified by hinge moment were formed, 


r ^pos“ V * IM 8sTALL^ 

^neg = " 

where the hinge moment, HM§, is positive for an external load 
resisting positive actuator motion, and HMg^ ^ represents 

the maximum load the actuators can overcome. Note that an 
aiding load will produce a rate limit larger than the no-load 
rate limit, rl^. An aiding load could occur for a leading edge 

surface. The positive and negative limits were imposed on the 


linear rate x 0[in 

as follows. 




r ! 

r *pos 

xo n„ > rlpos 

xo = 

< 

x o„„ 

otherwise 



r I 

l ' 4 neg 

ifx °/jn < rl ncg 


The state xq was then used as a command to the second-order 
portion of the actuator transfer function 

5 = (x 0 - 8) - 2 £$(0§ 5 

Once the position, velocity, and acceleration of the control 


Turbulence Model 

A turbulence model based on the Drydcn atmospheric 
turbulence model l 9 J was used. A break frequency of 
17.23 Hz for the turbulence transfer function was used to 
approximate the expected wind-tunnel turbulence. Resonance 
peaks at 10 Hz were observed in tunnel data from prior entries 
and 10-11 Hz was the predicted flutter frequency. A break 
frequency of 17.23 Hz produces a peak magnitude at 10 Hz in 
the Drydcn transfer function. The state equations used for 
each symmetry are, 

'7 3/2 1/2 

x g = -2(O g Xg - 0) g Xg + <Jg(0g' T u D 


where, \) is a Gaussian random process that is sampled and 
held every T\> seconds, (Og is the break frequency in 
radians/sccond, and Og is the RMS intensity. In both the 
batch and hot-bench simulations, T v is the integration step 

size. The factor of T^ applied to the input V arises from the 
inherently band-limited nature of a digital simulation. The 
random process that produces o can be interpreted as white 
noise being passed through a perfect (1/2T V ) Hz band-pass 
Filter. The output equations are, 


* ^ . 
5g = x g + -« g 


’fi 


. ; • V3 

a nd ^g = *g+— 

UJo 
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Prior to the November 1989 entry* an overall intensity 
level of 12 in/sec was estimated for the tunnel turbulence. It 
was decided to allocate the turbulence between the symmetric 
and anti -symmetric models based on an 85/15 percent 
distribution. The following symmetric and anti-symmetric 
intensities resulted 

°gsy = *°- 2in / sec °g a s = 1 * 8in / scc 

When the model is in the tunnel with a flutter-suppression 
control law engaged and subjected only to natural turbulence 
as excitation, control surface activity results. Analysis of data 
generated in the November 1989 tunnel test revealed that the 
RMS control surface activity predicted by simulation was 
much higher than actually experienced in the tunnel. Three 
different flutter control laws were tested in the November 
1989 tunnel test and all generally resulted in the same 
observed RMS levels and the same degree of ovcrprediciion 
by the simulation. By making the turbulence intensities (Gg) 
functions of dynamic pressure, simulation-predicted RMS 
levels can be brought into agreement with observed data. An 
example of this process is shown in figure 5 for the flutter 
suppression control law that achieved the highest dynamic 
pressure. The intensities required to bring batch-simulation- 
gencrated RMS results for both commanded control positions 
and control rates into agreement with experimental data arc 
plotted in figure 5 as functions of dynamic pressure. For each 
symmetry, one intensity function corrects the RMS rate 
predictions and the other corrects the RMS deflection 
predictions. Since only one intensity exists per symmetry, 
either the RMS rate results or the RMS position results, but 
not both, can be matched. The solid lines (one per symmetry) 
on figure 5 represent a fitted line biased upwards. The 
upward bias favors overprcdiction, which is conservative and 
preferable, if small. The solid lines on figure 5 are given by 
the following functions, 


°gsy 


= 0.4 + 



a Sas 


1.6 + 2.4 



where q has units of psf. These estimates are configuration 
dependent and should not be be regarded as a final 
characterization of turbulence in the TDT. 


Aeroclastic Equations 

The aeroclastic equations in a frequency domain format 
are given by HU 0,11] 


-m 


[ M f ] [M fc ] ‘ 
[M*] [ M 0 -] 


+ jeo 


[ G‘.] 0 

0 0 


[•K f ] 

0 


0 

0 




[Q ff (jco)l [Q fc (jco>] 1 M 
[Q cf Cjo»l fQ cc (jco)] J [Sj 


and apply to either symmetric or anti-symmetric motion. The 
flexible modes are augmented with control modes that 
represent control surface deflection. The control modes have 
zero stiffness. A low frequency subset of the elastic modes of 
free vibration from a large-order structural model are typically 
used in an aeroclastic formulation. For the AFW simulations, 
the number of retained flexible modes per symmetry was 
always between 7 and 10, inclusive. The flexible modes are 
orthogonal to each other so the mass and stiffness matrices are 
diagonal. Modal damping of 0.03 is assumed for each mode. 

The in-vacuo equations were augmented with generalized 
aerodynamic force coefficient matrices, [Q(j<o)L that are 
functions of frequency. The functions, [Q(jco)], can be 
approximated^ 0 ' 11 ^ by matrix expressions that are rational in 
the Laplace variable, s. The "least square"! 10 ’ 11 } form of 
approximation is given by 


[Q(s)l = [A 0 ] + fAjlxs + [A 2 ](xs)2 



( 2 ) 


where x = (c/2V) and is 1,2,3 or 4, depending on the 
order of the fit. Equations (1) and (2) can be combined to 
produce the time domain aeroclastic equations (3) in table 1, 
wherein the second-order, in-vacuo equations were 
augmented with unsteady aerodynamic "lag" states arising 
from the denominator term being summed in equation (2). 
Control surface positions, rates, and accelerations along with 
turbulence arc treated as external inputs. The vectors [r \ ) 

and | x a IT1 | arc both n f xl, where is the number of retained 

clastic modes. Equations (3) are used to solve for the clastic 
mode accelerations, {ti } , which can be integrated to find the 
rales and displacements. Equations (4) in table 1 are used to 
calculate actuator hinge moments. A positive hinge moment 
in this case resists positive actuator motion. The derivative 
calculations indicated in equations (3) and (4) were performed 
for each symmetry in the simulations. The symmetric and 
anli-symmctric components of the final accelerometer outputs 
were resolved into right and left components before output. 

Anti-Aliasing Filters 

For all the simulations, the dynamics of the anti-aliasing 
filters on the 40 primary outputs were simulated. For the 
tunnel test, both single-pole Fillers with a break frequency of 
25 Hz and fourth-order Butterworth Filters with break 
frequencies of 100 Hz, had been assembled and were 
available. Only the single-pole filters were used in the 
November 1989 tunnel test. The single pole anti-aliasing 
filters are given by. 



0) 


x(s) = 1 

U ( S ) (s/03 aa )+ 1 


where C0aa is ihe break frequency in radians/second. 


These equations consist of the standard in-vacuo second-order 
matrix structural equations (mass, damping, and stiffness) 
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Table 1 Aeroelastic and Hinge Moment Equations 



Pre-Test and Post-Test Models 

When the initial batch simulation was being developed, it 
was assumed the test would be conducted in Freon in the 0.8 
to 0.9 Mach range, and a real possibility existed of holding 
Mach constant during a test run and bleeding in Freon to raise 
the density of the test medium. It was also expected that some 
of the control laws would be scheduled on dynamic pressure 
in which case both the batch and hot-bench simulations would 
need to be able to vary dynamic pressure during a run to 


Table 2 Pre-Test and Post-Test Simulation Math Models 


State Categories 

Pre-test (1-lag) 

Post-test (4-lag) 

Symmetric flexible mode positions and velocities 

16 


20 


Sym aero lag states associated w/ the flexible modes 

8 


40 


Sym aero lag states associated w/ the control modes 

4 


16 


Symmetric turbulence states 

2 


2 


Anti-sym flexible mode positions and velocities 

14 


20 


Anti-sym aero lag states associated w/ the flexible modes 

7 


40 


Anti-sym aero lag states associated w/ the control modes 

4 


16 


Anti-sym turbulence states 

2 


2 


Total aeroelastic states 


57 


156 

Actuator states, 3 per actuator, 8 actuators 

24 


24 


Anti-aliasing Filters on 40 channels 

40 


40 


Total states 


121 


220 


effectively lest the control laws. The aerodynamic data arrays 
in the simulation are strong functions of Mach in the 0.8 to 
0.9 range. To interpolate each element of the aerodynamic 
data arrays according to Mach would require excessive CPU 
time. Therefore, the approach used in the both the pre-test 
batch and hol-bcnch simulations was to leave Mach and 
airspeed fixed and to vary the density to achieve a change in 
dynamic pressure. 
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Some months prior to the 1989 tunnel entry, the use of 
Freon as a test medium was forbidden, and it was determined 
that the test would be conducted with air as the test medium in 
the 0,2 to 0.5 Mach range. In this Mach range the 
aerodynamic matrices are virtually constant with respect to 
Mach. Furthermore, in the actual wind-tunnel test, the tunnel 
was not evacuated to any degree. Dynamic pressure, Mach, 
and airspeed were all changing as the fan speed was gradually 
increased. In the post-test implementation, density was left 
fixed, while airspeed was varied to replicate a given TDT 
dynamic pressure. 

As seen in table 2, the number of states in the post-test 
simulation math model is almost twice the number required in 
the pre-test math model. The post-test model uses 20 clastic 
modes instead of 15. For the post-test model, a 4-Iag least- 
square aerodynamic fit is employed (n/ a ^=4) instead of the 
l-lag fit used in the pre-test model. As seen in equations (3) 
and (4) and table 2, die choice of n{ ag has a dramatic impact 

on the number of states. 

Batch Simulation Implementation 

The structure of the batch simulation is identical to the 
simulation math model as described by differential equations. 
The state derivatives are all calculated explicitly from the 
states (outputs of integrators). The state derivatives arc 
collected in a vector and integrated with a Runge-Kutta 
second-order predictor -corrector method. The integration 
step used in the batch simulation is 1/2000 seconds. As 
indicated in figure 6, an integration step of 1/1600 seconds 
results in a small change in predicted response with 
significant degradation occurring for larger steps. In addition 
to the actuator, turbulence, aeroelastic, and filter dynamics, 
the batch simulation also simulates the digital controller. The 
effects of computational delay and quantization arc modeled. 

Pre-Test Hot-Bench Implementation 

The pre-test hot-bench simulation was implemented in a 
manner very similar to the batch simulation. State vectors 
(and associated state derivative vectors) that included all but 



the anti-aliasing filter states were formed. The vector of state 
derivatives was integrated numerically with no assumptions 
of linearity. Once current-time flexible mode accelerations 
were calculated from the current-time positions and velocities, 
an Adams-Bashforth second-order (AB2) predictor method 
was used to predict the velocities at the next time step. These 
predicted velocities were then used in a trapezoidal integration 
scheme to generate predicted flexible mode positions. To wit: 

^l(t+h) = n(t) + ^ ( 3 rj(t) - ri(t-h) ) 

Ti(t+h) = ri(t) + |( n(t+h) + Ti(t)) 

Using this modified AB2-based method, accuracy comparable 
to the Rungc-Kulta second-order predictor-corrector formula 
(RK2) used in the batch simulation was achieved with a 
single-pass formula. Note that the RK2 formulation requires 
two derivative evaluations per time step, whereas AB2 is a 
single-pass formula which gives up some gain accuracy for 
phase accuracy and is typically used in real-time applications. 
The integration step of 1/2000 seconds used in the pre-test 
hot-bcnch simulation was small enough that the batch and 
hot-bcnch simulation results compared favorably. 

The anti-aliasing filters were handled separately from the 
aeroelastic and actuator states. A scalar form of the state 
transition equations for a constant input over the interval was 
used. For the single-pole anti-aliasing filters given by, 

x(s) _ 1 

u(s) (s/CDaaH 1 

the constant input state transition equation is 

x(t+h) = c" ((0aah) x(t) + (l-e' (0)aah) )u(t) 

As seen in equations (3), calculating the accelerations of 
the clastic modes requires solving a matrix equation involving 
the mass matrix. If the mass matrix is constant, an inversion 
can be performed prior to a run (in a "reset" mode) and the 



o 0.1 0.2 

Time, seconds 


Figure 6.- Effect of step sizes on batch simulation response. 
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results simply stored for use as the integration proceeds. The 
mass matrix for the flexible modes, T M -]» (see equation (1)) 

is both constant and diagonal, and its inversion is trivial. The 
s-plane formulation used for the unsteady aerodynamic states, 
however, augments the mass matrix with a fully populated 

ff 

matrix, [A 2 L that is a function of Mach number. The total 
mass matrix becomes 

[Ml = [M f ] + p (j) [A*] 



(5) 


= ([!] + [A])['M f ] 


If the density of the test medium, p, is to be varied from the 
simulation control console without losing time 
synchronization, then the mass matrix equation must be 
solved at each derivative evaluation. Because the elements of 


[A] in equation (5) were much smaller than unity, the 
following approximation was successfully used for the 
inverse of the mass matrix, 

*2 ~ - 'N 


[M]l =['M f J-! 


[I] - p( j) [A f 2 f ]['M f .]-1 


When checked against the exact answer for maximum 
anticipated values for p, induced errors for the pre-test model 
were less than the precision available in the analog/digital 
converters used in the simulation loop. 

Prior to the November 1989 test, an integration step size 
of 1/2000 seconds and a compute frame of 1/80 seconds were 
required for the hot-bench simulation. A compute frame of 
1/80 seconds was achieved only if the AFW simulation was 
the only job on a Cyber 175. Thus, at best, the hot-bench 
simulation ran 25 (2000/80) times slower than real time, i.e., 
at a time scale of 1:25. A 1:25 time scale, however, proved to 
be burdensome while testing the controller performance 
evaluation mode. Sine sweeps taking 20 seconds in the 
tunnel would take over 8 minutes in the hot-bench lab. 
Modifications to the hot-bench simulation implementation, 
discussed in the next section, allowed the simulation to use a 
1/400 second integration step, allowing a 5-fold improvement 
in time scale. 


The post-test implementation method, wherein the hot- 
bench simulation is updated by data extracted from the batch 
simulation, is depicted schematically in figure 7. The second- 
order part of the third-order actuator models can be lumped 
with the remaining linear dynamics. The box in figure 7, 
labeled "Acroclaslic, Hinge Moment, and Load Equations," 
represents the slate equations of table 1 together with 
algebraic output equations to estimate the required 
accelerometer outputs, strain gage outputs, and pressure 
transducer outputs from the stales. Together with the 16 
stales associated with the second-order part of 8 actuator 
transfer functions and N (57 pre-test, 1 56 post-test) states 
from the aeroelaslic model, a coupled linear system of N+16 
states, 10 inputs (8 actuator and 2 noise), and 40 outputs can 
be extracted from the "linear" portion of the batch simulation. 
To implement the pre-test math model using the post-test state 
transition method required no model reduction on the 
extracted model, i.e., the simulation calculations could be 
completed in the same 1/80 second real-time clock frame used 
by the pre-test simulation, resulting in a time scale of 1:5. To 
maintain a 1 :5 time scale ratio while implementing the large 
post-test math model will require model reduction techniques 
to be applied to the extracted model. Reduction methods 
utilizing internal balancing techniques are being investigated. 
Once the model has been reduced, the state transition model 
based on an integration step of 1/400 seconds is calculated. 

The nonlinear portion involves only eight states, one from 
each actuator. Each state is integrated numerically with an 
integration step of 1/1600 seconds. Four integration steps are 
made to predict the value of the input to the coupled linear 
system at time (k+l)h where h = 1/400 seconds. Since input 
to the coupled linear system at time (k+l)h is now available, a 
trapezoidal state transition scheme can be employed. Let 
(u k ) denote the quantity (u(kh)). Given the linear dynamic 

system 


(x) = [A] (x) + [B]{u) 
if the ramp input signal, 

(u(t)} = K) + (t-kh) frfriHH kl 

is defined over the interval 

kh <t<(k+l)h 

then the following exact solution for (x) at time t = (k+l)h 
exists 


Post-Test Hot-Bench Implementation 

The post-test hot-bench simulation implementation was 
driven by the need to reduce the time-scale factor to 
something closer to real time together with the need to 
accommodate the larger post-test models. If the hot-bench 
simulation is restricted to a fixed tunnel operating point for a 
given run, i.e., density, Mach, and airspeed are held fixed, 
then once the rate limiting is performed on the actuator 
transfer functions, the remaining dynamics in the simulation 
are linear. By utilizing a state transition method of 
discretization on these dynamics, the hot-bench integration 
step has been increased by a factor of 5, from 1/2000 to 1/400 
seconds. 


{x k+ ,} = [F]{x k ] + [G 0 ](u k ] +[G,]{u k+1 ) 

where, 

m = * lA]h 

[C 0 1 = (c rAlh |Ahl' 1 - [Ah] 1 - e [Alh )[-A]' , [B] (6) 

[G,] = ([I] - e IAlh [Ah]' 1 + [AhrOl-Ar’PJ (7) 
Note that 

[G] = [G 0 ] + [G,] = ([I] - e [A]h ) [-A] _1 [B] 
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Batch Simulation 



Figure 7.- Data flow from the baleh to the hot-bench simulation. 


which corresponds to the result one expects to see for [G] if 

u(t) is assumed to be constant over the interval ([A]h)c^^ = e^^[A]h 


kh < t<(k+l)h 

Clearly, the direct evaluation of [Gq] and [Gj] using 
equations (6) and (7) will not work if [A] is singular, as 
would occur if [A] included rigid body modes with zero 
eigenvalues. However, using the Taylor series expansion 

e M h = [I] + [A]h + A<[A]h)2 + _ 

oo 

= X 5! <[Alh,P 

p=0 


and recognizing that 


the equations for [Gq] and [Gj] can be put into a form that 
can be calculated if [A] has zero eigenvalues. Thus, 


( 


[G 0 ] = 


^ ^j(-l) P ([A]h)P' 2 
V P=2 


h e [A]h [B] 


[Gil = 


oc 

'Lh. 


( 8 ) 


«A]h) 


P-2 


h [B] 


V P=2 
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The matrices [Gq] and [Gi] can be calculated by summing the 
above series until the next term is under some tolerance. 

When applied to the pre-test model, procedures to sum the 
series defined by equations (8) converged without difficulty. 

The anti-aliasing filters are applied individually to each 
output signal, which results in a diagonal system. The anti- 
aliasing filters are therefore not lumped with the coupled 
linear system to avoid making full matrix operations. The 
anti-aliasing filter dynamics are digitized in a sequential 
manner utilizing a scalar form of the trapezoidal state 
transition method described above. For single pole anti- 
aliasing filters given by 

x(s) _ 1 

u(s) (s/cO+1 

the state transition equations are 

x(t+h) = e' (C0aah) x(t) + gou(t) + g,u(t+h) 

where 

go = -e' (<0 ^ h Wh)-i + (oWO > - c' (w - h) (9) 

g, = 1 + e^W)- 1 - (co^h)-! (10) 

Note that the term (-[A]' 1 [B]) in (6) and (7) becomes unity 
in equations (9) and (10). 

Verification 

The hot -bench simulation was verified by comparing time 
history results with the batch simulation. Trajectories were 
found to match within the width of plotted lines. The open- 
loop plant dynamics of the batch simulation were verified by 
comparison with IS AC-generated linear models. For each 
symmetry, a linear model was extracted from the batch 
simulation using finite differencing. A batch-generated 
symmetric model and an ISAC-generated symmetric model 
had different numbers of states because of the way the 
actuators were handled. The batch-derived model had 
dynamics for eight right and left control surfaces and the 
ISAC model had dynamics for only four symmetric control 
deflections. The corresponding batch and ISAC linear 
models were compared by overlaying the gain and phase 
frequency response of each input/output pair. The agreement 
between ISAC-generated and batch-simulation-dcrivcd 
frequency responses was excellent. 

Concluding Remarks 

Two simulations, one batch and one real-time, of an 
aeroelastically scaled wind-tunnel model were described. The 
batch simulation was used to generate and verify the real-time 
simulation and to test candidate control laws prior to 
implementation on the control computer. The real-time 
simulation supported hot-bench testing of a digital controller 
developed to actively control the elastic deformation of the 
wind tunnel model. Time scaling required for hot-bcnch 
testing was discussed. Substantial improvement in the time 
scale ratio was achieved by application of state transition 
methods. 
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Appendix A - Input/Ouput List for AFW Wind-Tunnel Model and Simulations 


Tabic A.l Simulation and Wind-Tunnel Model Inputs and Scalin 


No. 

Signal 

Volts/Um 

1 Units 

1 

Left LEO actuator command, + leading edge down 

0.375 

degrees streamwise 

2 

Left LEI actuator command, + leading edge down 

0.375 

degrees streamwise 

3 

Right LEI actuator command, + leading edge down 

0.375 

degrees streamwise 

4 

Right LEO actuator command, + leading edge down 

0.375 

degrees streamwise 

5 

Left TEO actuator command, + trailing edge down 

0.375 

degrees streamwise 

6 

Left TEI actuator command, + trailing edge down 

0.375 

degrees streamwise 

7 

Right TEI actuator command, + trailing edge down 

0.375 

degrees streamwise ? 

8 

Right TEO actuator command, + trailing edge down 

0.375 

degrees streamwise 


Table A. 2 Simulation and Wind-Tunnel Model Outputs and Scaling 


No. 

Signal 

Volts/Unit 

Units 

1 

Left LEO actuator position, RVDT, + leading edge down 

0.375 

degrees streamwise 

2 

Left LEI actuator position, RVDT, + leading edge down 

0.375 

degrees streamwise 

3 

Right LEI actuator position, RVDT, + leading edge down 

0.375 

degrees streamwise 

4 

Right LEO actuator position, RVDT, + leading edge down 

0.375 

degrees streamwise 

5 

Left TEO actuator position, RVDT, + trailing edge down 

0.375 

degrees streamwise 

6 

Left TEI actuator position, RVDT, + trailing edge down 

0.375 

degrees streamwise 

7 

Right TEI actuator position, RVDT, + trailing edge down 

0.375 

degrees streamwise 

8 

Right TEO actuator position, RVDT, + trailing edge down 

0.375 

degrees streamwise 

9 

Model pitch actuator position, RVDT, 4* nose up 

0.375 

degrees 

10 

Model roll position, + right wing down 

0.0555 

degrees 

11 

Left LEO co-located accelerometer, + up 

0.5 

g 

12 

Right LEO co-located accelerometer, + up 

0.5 

g 

13 

Left TEO co-located accelerometer, + up 

0.5 

g 

14 

Left TEI co-located accelerometer, + up 

0.5 

g 

15 

Right TEI co-located accelerometer, + up 

0.5 

g 

16 

Right TEO co-located accelerometer, + up 

0.5 

g 

17 

Left wingtip accelerometer, + up 

0.5 

g 

18 

Right wingtip accelerometer, + up 

0.5 

g 

19 

Left store-mounted accelerometer, + up 

0.5 

g 

20 

Right store-mounted accelerometer, + up 

0.5 

g 

21 

Fuselage accelerometer # 1, + up 

1.0 

g 

22 

Fuselage accelerometer # 2, + up 

1.0 

g 

23 

Fuselage accelerometer # 3, + up 

1.0 

g 

24 

Roll rate, + right wing down 

0.0224 

deg/sec 

25 

Left outboard bending moment, + tip up 

0.00244 

in-lb 

26 

Left inboard bending moment, + tip up 

0.000477 

in-lb 

27 

Right inboard bending moment, + tip up 

0.000553 

in-lb 

28 

Right outboard bending moment, + tip up 

0.002820 

in-lb 

29 

Left outboard torsion moment, + leading edge up 

0.00611 

in-lb 

30 

Left inboard torsion moment, + leading edge up 

0.000112 

in-lb 

31 

Right inboard torsion moment, + leading edge up 

0.000106 

in-lb 

32 

Right outboard torsion moment, + leading edge up 

0.00702 

in-lb 

33 

Left LEO actuator hinge moment, + leading edge up 

0.014760 

in-lb 

34 

Left LEI actuator hinge moment, + leading edge up 

0.014144 

in-lb 

35 

Right LEI actuator hinge moment, + leading edge up 

0.014155 

in-lb 

36 

Right LEO actuator hinge moment, + leading edge up 

0.020503 

in-lb 

37 

Left TEO actuator hinge moment, + trailing edge up 

0.026917 

in-lb 

38 

Left TEI actuator hinge moment, + trailing edge up 

0.014592 

in-lb 

39 

Right TEI actuator hinge moment, + trailing edge up 

0.013616 

in-lb 

40 

Right TEO actuator hinge moment, + trailing edge up 

0.028341 

in-lb 

41 

Wind Tunnel Mach number 

10.0 

Mach number 

42 

Wind Tunnel Dynamic Pressure 

.006945 

PSI 
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